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ABSTRACT 

Collisional growth of submicron-sized dust grains into macroscopic aggregates is the first step of planet for- 
mation in protoplanetary disks. These grains are expected to carry nonzero negative charges in the weakly 
ionized disks, but its effect on their collisional growth has not been fully understood so far. In this paper, we 
investigate how the charging affects the evolution of the dust size distribution properly taking into account 
the charging mechanism in a weakly ionized gas as well as porosity evolution through low-energy collisions. 
To clarify the role of the size distribution, we divide our analysis into two steps. First, we analyze the colli- 
sional growth of charged aggregates assuming a monodisperse (i.e., narrow) size distribution. We show that the 
monodisperse growth stalls due to the electrostatic repulsion when a certain condition is met, as is already ex- 
pected in the previous work. Second, we numerically simulate dust coagulation using Smoluchowski's method 
to see how the outcome changes when the size distribution is allowed to freely evolve. We find that, under 
certain conditions, the dust undergoes bimodal growth where only a limited number of aggregates continue to 
grow carrying the major part of the dust mass in the system. This occurs because remaining small aggregates 
efficiently sweep up free electrons to prevent the larger aggregates from being strongly charged. We obtain a 
set of simple criteria that allows us to predict how the size distribution evolves for a given condition. In Paper 
II, we apply these criteria to dust growth in protoplanetary disks. 

Subject headings: dust, extinction — planetary systems: formation — planetary systems: protoplanetary disks 



1. INTRODUCTION 

The standard core-accretion scen ario for planet formation 
dMizunolll980t iPollack et al.ll 19961 ) is based on the so-called 
planetesimal hypothesis. This hypothesis assumes that solid 
bodies of size larger than kilometers (called "planetesimals") 
form in a protoplanetary disk prior to planet formation. How- 
ever, the typical size of solid particles in interstellar spac e 
is as small as a micron or even smaller (Mathis et al.ll 1977b . 
It is still open how the submicron-sized grains evolved into 
kilometer-sized planetesimals. 

The simplest picture for dust evolution towards planetesi- 
mals can be summarized into the following steps. (1) Initially, 
submicron-sized particles coagulate into larger but highly 
porous, fractal aggregates through low-velocity collisions 
driven by Brownian motion and differential settl ing towards 
the midplane of th e disk (IWurm & Bluml 119981: iBlum et all 
Il998t iKempf et all 119991) . (2) As the aggregates grow to 
"macroscopic" (mm to cm) sizes, the collisional energy be- 
comes high enough to cause the compaction of the aggregate s 
(Blum 2004| ISuvamaetal] 120081: iPaszun & Dominiki 120091) . 
(3)The compaction cause the increase in the stopping times 
of the aggregates, allowing them to concentrate in the mid - 
plane of the disk dSafronovll 19691 IGoldreich & Ward fl973l 
the center o f vortices dBarge & Sommerial 11995b . or turbu- 
lent eddies ( Johans en et all 120071) . (4) Planetesimals may 
form within such dense regions through gravitational insta- 
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bility (ISafronovll 19691 IGoldreich & Wardl [19731) or through 
further collisional gr owth jWeidenschilling & Cuzzil 1 1993b 
IWeidenschilhnglll995l) . 

However, there is great uncertainty on how large dust 
aggregates can grow through mutual collisions (see, e.g., 
IBlum & Wunnll200llGuttler et alj|2010l) . As the collisional 
compaction proceeds, the aggregates decouple from the am- 
bient gas and obtain higher and higher re lative velocities 
driven by radial drift (|Weidenschillinglll977[) and gas turbu- 
lence dVolk et alj|198u) . The collision velocity can exceed 
10 ms" 1 even without turbulence, but it is uncertain whether 
such high-speed collision s lead to the sticking or fragmen- 
tation of the aggregates dBlum & Wurml l2008: Wad a et alj 
120091 iTeiser & Wurml 120091: IGuttler et all 120101) . In addi- 
tion, collisional co mpaction itself can cause the reduction o f 
sticking efficiency dBlum & Wurm 12001 iGiitfler et alj20Tob . 
This may termina tes the collisiona l growth before the frag- 
mentation occurs dZsom et al.ll2010l) . 

By contrast, it is generally believed that dust coagulation 
proceeds rapidly until the aggregates grow beyond the initial, 
fractal growth stage since the collision v elocity is too low to 
cause the reduction of sticking efficiency (jD ominik & Tielensl 
[19971 IBlum & Wurml 12001 IGuttler et alj|2010t) . However, 
one of the authors has recently pointed out that electric charg- 
ing of aggregates cou ld halt dust gro wth before the aggre- 
gates leave this stage (Okuzumi 2009, hereafter O09). Pro- 
toplanetary disks are expected to be weakly ionized by a 
various kinds of high-energy sources, such as cosmic rays 
( Umebavashi & Nakano 1981) and X-rays from the central 
star ([Glassgold et al. 1997). In such an ionized environment, 
dust particles charge up by capturi ng ions and electrons, a s 
is well known in plasma physics (Shukla & Mamunl 120021) . 
In equilibrium, dust particles acquire nonzero negative net 
charges because electrons have higher thermal velocities than 
ions. This "asymmetric" charging causes a repulsive force 
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between colliding aggregates, but this effect has bee n igno red 
in previous studies on protoplanetary dust growth. IO09I has 
found that the dust charge in a weakly ionized disk can be con- 
siderably smaller than in a fully ionized plasma but can never- 
theless inhibit dust coagulation in a wide region of the disk. It 
is also found that the electrostatic barrier becomes significant 
when the dust grows into fractal aggregates, i.e., much earlier 
than the growth barriers mentioned above emerge. Thus, the 
dust charging can greatly modify the current picture of dust 
evolution towards planetesimals. 

The analysis of the electrostatic barrier by O09 is based on 
the assumption that dust aggregates obey a narrow size distri- 
bution. In reality, however, size distribution is determined as a 
result of the coagulation process, and it has been unclear how 
the distribution evolves when the dust charging is present. The 
purpose of this study is to clarify how the size distribution of 
dust aggregates evolves when the aggregates are charged in a 
weakly ionized gas. 

According to O09, the effect of dust charging can be- 
come already significant before the collisional compaction 
of aggregates becomes effective. In this stage, dust aggre- 
gates are expected to have lower and lower internal density 
(i.e., higher and higher porosity) as they grow, as is sug- 
gested by laboratory experiments and j V-body simulation s 
dWurm & Blum|[T99l IBlum et alJ[l998t IKempf et al]fl999l) . 
This porosity evolution has been igno red in most theoreti- 
cal studies on dust coagulation (e.g., iN akagawa et al.l [198 lb 
iTanaka et al] 120051 [Brauer et al. 2008), in which aggregates 
are simplified as compact spheres. However, when analyz- 
ing the electrostatic barrier, the porosity evolution must be 
accurately taken into account; in fact, as we will see later, 
the ignorance of the porosity evolution leads to considerable 
underestimation of the electrostatic barrier, because compact 
spheres are generally less coupled to the ambient gas and 
hence have higher collision energies than porous aggregates. 
In this study, we use the fractal dus t model recently proposed 
bv lOkuzumi et al.l ( 120091 hereafter lOTS09l) . Classically, frac- 
tal dust growth has been only modeled with either of its two 
extreme limits, namely, ballistic cluster-clus ter and particle- 
cluster aggregation (BCCA an d BPCA; e.g..rOssenkop fll993l: 
Dullemond & Dominik 2005). To fill the gap between the two 
limits, IOTS09I introduced a new aggregation model (called 
the quasi-BCCA model) in whic h aggregates grow through 
unequal-sized collisions. IOTS 09 found from iV-body simu- 
lations that the resultant aggregates tend to have a fractal di- 
mension D close to 2 even if the size ratio deviates from unity. 
This explains why fractal aggregates with D ~ 2 are univer- 
sally observed in various low-velocity coagulation processe s 
(Wur m &Bluml[T998t IBlum et alJ[l998t IKempf et al.l[l999l) . 
IOTS09I summarized the results of their A^-body simulations 
into a simple analytic formula giving the increase in the poros- 
ity (volume) for general hit-and-stick collisions. This formula 
together with the Smoluch owski equation extended for porous 
dust coagulation (OTS09) enables us to follow the evolution 
of size distribution and porosity consistently with dust charg- 
ing. 

As we will see later, our problem involves many model pa- 
rameters, such as the initial grain size and the gas ionization 
rate. To fully understand the dependence of the results on 
these parameters, we do not assume any protoplanetary disk 
model but seek to find general criteria determining the out- 
come of dust evolution. This approach allows us to investigate 
the effect of the electrostatic barrier with any protoplanetary 
disk models. Application of the growth criteria to particular 




FIG. 1. — Projection of a numerically created, three-dimensional porous 
aggregate consisting of ss 1000 monomers. The large open circle shows the 
characteristic radius a (for its definition, see Section 2.3.1), while the gray 
disk inside the circle shows the projected area A averaged over various pro- 
jection angles. Note that A is not necessarily equa l to Tea 2 , especially when 
the aggregate is highly porous (see also Figure 4 of OTS09). 

disk models will be done in lPaper fj dOkuzumi et al.ll20TTl) . 

This paper is structured as follows. In Section 2, we de- 
scribe the dust growth model used in this study. In Section 3, 
we examine the case of monodisperse growth in which all 
the aggregates grow into equal-sized ones. The monodis- 
perse model allows us to introduce several important quan- 
tities governing the outcome of the growth. We analytically 
derive a criterion in which the "freezeout" of monodisperse 
growth occurs. In Section 4, we present numerical simula- 
tions including the evolution of the size distribution to show 
how the outcome of the growth differs from the prediction of 
the monodisperse theory. We discuss the validity of our dust 
growth model in Section 5. A summary of this paper is pre- 
sented in Section 6. 

2. DUST GROWTH MODEL 

In this section, we describe the dust growth model consid- 
ered in this study. 

We consider collisional growth of dust starting from an en- 
semble of equal-sized spherical grains ("monomers"). Each 
aggregate is characterized by its mass, radius, projected area, 
and charge. For simplicity, we assume "local" growth, i.e., we 
neglect global transport of dust within a disk. 

We focus on the first stage of dust evolution in protoplan- 
etary disks and assume that aggregates grow through "hit- 
and-stick" collisions, i.e., collisions with perfect sticking ef- 
ficiency and no co mpaction. It is known th eoretically (e.g. 
Kemp f et al.l 119991) and experimentally (e.g. Wurm & Bluml 
1998) that hit-and-stick collisions lead to highly porous aggre- 
gates. To take into account the po rosity ev olution, we adopt 
the fractal dust model proposed by OTS09. This model char- 
acterizes each a ggrega te with its mass M and "characteristic 
radius" a (see OTS09 and Section 2.3 for the definition of 
the characteristic radius), and treat the two quantities as in- 
dependent parameters. Another important parameter is the 
projected area A This determines ho w the ag gregates are fric- 
tionally coupled to the gas. In the OTS09 model, A is not 
treated as an independent parameter but is given as a function 
of M and a. Note that A is not generally equal to a naive "cross 



ELECTROSTATIC BARRIER AGAINST DUST GROWTH. I 



3 



section" na 2 , especially when the aggregates is highly porous 
(FigureQ"! see also Figure 4 of OTS09). Distinction between 
A and ira 2 allows us to avoid overestimation of the gas drag 
force to dust aggregates. In Section 2.3, We will describe the 
porosity model in more detail. 

The collision probability between two aggregates 1 and 2 
is proportional to their relative s peed Au times the collisio nal 
cross section <7 co u given by (e.g., Landau & Lifshitz 1976) 



ion-electron plasma (IEP) 



I n(a\ + C12) I 1 ■ 

Ccoll = \ 

0, 



Eel 
^kin 



^kin ^ ^el 
^icin ^ ^el j 



(1) 



where E^ n = M M ( A«) 2 /2 is the kinetic energy associated with 
the relative motion, M M = M l M 2 /(M l +M 2 ) is the reduced 
mass, and E e \ = Q\Q2/{a\+ai) is the energy needed for the 
aggregates to collide with each other. In this paper, E e \ is 
called "the electrostatic energy" for colliding aggregates. Be- 
low, we describe how to determine Q and Au. 

2.1. Charging 

We adopt the dust charging model developed by IO091 In 
this model, dust aggregates are surrounded by a weakly ion- 
ized gas and charge up by capturing free electrons and ions. 
These ionized particles are created by the nonthermal ioniza- 
tion of the neutral gas and are removed from the gas phase 
through the adsorption to the dust as well as the gas-phase re- 
combination. The dust charge Q and the number densities of 
ions and electrons are thus determined by the balance among 
the ionization, recombination, and dust charging. In equilib- 
rium, the average charge (Q) a of aggregates with radius a is 
given by (see Equation (23) of lO09l) 



(<2} fl = -* 



aknT 



(2) 



where k% is the Boltzmann constant, T is the gas temperature, 
e is the elementary charge, and ^ is a dimensionless param- 
eter characterizing the charge state of the gas-dust mixture. 
O09 has analytically shown that the equilibrium conditions 
are reduced to a single equation for "J. When the adsorption 
to the dust dominates the removal of the ionized gas, the equa- 
tion for ^ is written as (see Equation (34) of O09) 



m,. , >l' „ 
-exptf-- = 0, 
nti <a 



(3) 



where m,-^ is the mass of ions (electrons), s;( e ) is their sticking 
probability onto a dust monomer, and 



6 = 



TTWj 



(4) 



is a dimensionless quantity depending on the total pro- 
jected area A tot = J A(M)n(M)dM and total radius C tot = 
J a(M)n{M)dM of aggregates, and the ionization rate £ and 
number density n g of neutral gas particles. Equation (01 orig- 
inates from the quasi-neutrality condition, en— en e + Q tot = 0, 
where n,- and n e are the number density of ions and electrons, 
and 2tot = J (Q)a(M)n{M)dM is the total charge carried by dust 
in a unit volumeQ Equation (f3]i cannot be used when the gas- 
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FIG. 2. — Schematic illustration of an ion-electron plasma (IEP: left) and 
an ion-dust plasma (IDP; right). In an IEP, the dominant carriers of negative 
charges are free electrons. In an IDP, by contrast, the dominant negative 
species is the charged dust. The absolute value of the dust surface potential, 
\il>\ =a|Q|, is generally smaller in IDPs than in IEPs. 



phase recombination dominates the removal of the ionized 
gas. In a typical protoplanetary disk, however, the gas-phase 
recombination can be safely neglected unless the dust-to-gas 
ratio is many orders of magnitude smaller than interstellar val- 
ues - 0.01 (IO09h . 

Physically, is related to the surface potential of aggre- 
gates. For an aggregate with charge Q and radius a, the sur- 
face potential ip is given by ip = Q/a. Equation (f2]i implies that 
^ = {^) a / {-k-eJ /e), namely, ^ is the surface potential aver- 
aged over aggregates of radius a and normalized by -k%T je. 
Note that (ijj) a is apparently independent of a, but is actually 
not because "J depends on the size distribution of aggregates 
through A tot and C tot . It should be also noted that the radius a 
can be interpreted as the electric capacitance C (i.e., Q = Cip). 
This is the reason why we have denoted the total radius as C tot . 

As sho wn in IO09L ^ asymptotically behaves as (see Sec- 
tion 2.3 of [O09) 



where is the solution to 



1 



■exp^oo = 0. 



(5) 



(6) 



Equation © is known as the equation for the equilibrium 
charge of a dust particle e mbedded in a fully ionized plasma 
dSpitzerlll94lt IShukla & Mam un 2002). Equation © sug- 
gests that the charge state of dust particles in a weakly ionized 
gas is characterized by two limiting cases. If 3> ^>oo, the 
total negative charge \Q to t\ carried by dust aggregates is neg- 
ligibly small compared to en e , and the quasi-neutrality con- 
dition approximately hold in the gas phase, i.e., n, w n e . If 
8 <C ^00, by contrast, most of the negative charge in the sys- 
tem is carried by aggregates, and the quasi-neutrality condi- 
tion approximately holds between ions and negatively charged 
dust. For this reason, O09 referred to the former phase as the 
ion-electron plasma (IEP), and to the latter as the ion-dust 
plasma (IDP). Figure [2] schematically shows the difference 
between the two plasma states. 

For given m, and Si/s e , Equation I© determines W as a func- 
tion of 9. In typical protoplanetary disks, the dominant ion 
species are molecular ions (e.g., HCO + ) or metal ions (e.g., 
Mg + ) depending on the abundance of metal atoms in the gas 
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FIG. 3. — Comparison between the numerical solutions to Equation {5J 
and the approximate formula (7)- The symbols indicate the numerical solu- 
tions for various values of s e , and the solid curves show the prediction from 
Equation Q. The ion mass is taken to be 24#ih for all the cases. The max- 
imum values are 3.78, 2.81, 1.96, and 1.10 for s e = 1.0, 0.3, 0.1, 0.03, 
respectively. 



phase dSano et al.ll2000t lllgner _& Nelsonll2006l). Although Sj 
is likely to be close to unity ( Um ebavashi & Nakanoi [19801; 
iDraine & Sutinll9"87l). S*. a t low temperatures is poorly under- 
stood. lUmebavashil d 19831) estimated s e using a semiclassical 
phonon theory to obtain 0.1 < s e < 1 for T < 100K. However, 
the uncertainty in s e does not strongly affect the evaluation of 
"J*. For example, assuming m, = 24m H (the mass of Mg + ) and 
Si = 1, is 3.78 for s e = 1, and is 1.96 even for s e = 0.1. 

Figure|3]illustrates the dependence of on 6 for fixed m,(= 
24m H ) and s,(= 1) with various s e (=l.0, 0.3, 0.1, 0.03). We 
find that W can be well approximated by 



1 + 



_e_ 



-1/0.8 



(7) 



In Figure [3] we compare Equation (|7]i with the numerical so- 
lutions to the original equation. The approximate formula 
recovers all the numerical solutions within an error of 20%. 
This means that ^/^oo is well approximated as a function of 
Q/^oo for this parameter range0 We use this fact in Section 
3. 

Up to here, we have considered only the mean value of the 
charge Q. In fact, there always exists a finite value of the 
charge dispersion (5Q 2 ) a , and moreover, the mean value (Q) a 

is not necessarily larger than {8Q 2 ) l J 2 (O09). Nevertheless, 
we will assume below that the dust charge Q is always equal 
to (Q)a- The validity of this assumption will be discussed in 
Section 5.2. 

2.2. Dust Dynamics 

As found from Equation ((T), the relative velocity between 
aggregates determines whether they can overcome the elec- 
trostatic barrier to collide. In this study, we model the motion 
of dust aggregates in the following way. We assume that the 
motion of each aggregate relative to the ambient gas consists 
of random Brownian motion and systematic drift due to spa- 
tially uniform acceleration (e.g., uniform gravity). Then, the 
probability density function P r (Au) for the relative velocity 

6 This is not true for more general cases. In fac t, Equations {3) and (6) can 
be combined into a single equation (Equation fl33t ). which cannot be reduced 
to an equation for ^/^oo depending only on O/'I'oo. 



Au = Ui -U2 between two aggregates 1 and 2 is given by 

MJAu-Au D ) 2 ^ 



P r (Au)dAu = ( 



M„ 



\2Trk B T 



,3/2 

J exp 



2k n T 



dAu, 



(8) 

where Aud is the difference of the drift velocities between 
the two aggregates. Here, we have assumed that the system- 
atic motion has no fluctuating component, that is, the velocity 
dispersion is thermal even when M^Au D 3> k B T. We will dis- 
cuss the effect of adopting a different velocity distribution in 
Section 5.3. 

We further assume that aggregates are frictionally coupled 
to the ambient gas, and give Aud as 



Au D = g(n-T 2 ), 



(9) 



where Tj(j = 1,2) is the stopping time of each aggregate and 
g is the uniform acceleration. In this study, we focus on small 
aggregates and give r according to Epstein's law, 



3M 



(10) 



where p g is the gas density and m g is the mass of the gas par- 
ticles. Epstein's law is valid when the size a of the aggregate 
is smaller the mean free path I of gas particles. 

In a protoplanetary disk, relative motion like Equation (0 
is driven by several processes. For example, the gravity of the 
central star causes acceleration g = Vt\z towards the midplane 
of the disk, where is the Kepler rotational frequency and 
z is the distance from the midplane. Another example is the 
acceleration driven by gas turbulence in the strong coupling 
limit. When both of two colliding aggregates are frictionally 
well coupled to the turbulent eddies of all scales, the rela- 
tive velocity between the aggregates is approximately given 
by Auo « (u v /t v y\T\— ti\, where u v and t v are the characteris- 
tic vel ocity and turnover time f or the smallest eddies, respec- 
tively dWeidenschillingl Il984t lOrmel & CuzzH 120071) . This 
means that turbulence behaves as an effective acceleration 
field of g « u v /trj for strongly coupled aggregates. 

As the collisional cross section <7 co u depends on the stochas- 
tic variable Au, it is useful to treat collision events statisti- 
cally. To do so, we introduce the collisional rate coefficient 



K 



P,(Au)cr co u|Au|c/Au. 



(11) 



With Equations ([TJ and ©, the integrati on can be a nalytically 
performed. Using g, Q 2 > 0, we have dShullll 1978b 



K=n(a\ +a 2 ) 



k K T 



2-kM^Sd 



y+ exp(-y^) - y_ exp(-y+) 



+■^-(1 -2y + y_){erf(y + )-erf(y_)} 



(12) 



where erflj) = (2/ ^/n) J? e\p(-z 2 )dz is the error function, and 
y + and y- are defined as 



with 




(a\+a2)k B T 



(13) 
(14) 
(15) 
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Note that £q and £ E are the relative kinetic energy associated 
with differential drift and the electrostatic energy normalized 
by k^T, respectively. 

Equation Sl2\ has the following simple asymptotic forms: 

!7r(ai +a 2 ) 2 Au B exp(-£ E ), £ D < 1, 
, A / £ E \ (16) 
7r(ai+a 2 ) Am d 11 -— I , £ D >1,£ £ , 

where Amb = {%k^T /-kM^) 1 ! 1 is the mean thermal speed 
between the colliding aggregates. The exponential factor 
exp(-£ e) originates from the high-energy tail of the Maxwell 
distribution. This factor guarantees K nonvanishing even for 
large £ E . 

2.3. Porosity Model 

As shown by IO091 the charging affects dust growth 
before the collisional compaction becomes effective. In 
this early stage, aggregates have a high ly porous structure 
dWurm & Blumlll998t iKempfet al.lll999l The porosity in- 
fluences their collisional growth through the collisional and 
aerodynamical cross sections. It also affects dust charging 
through the capacity (=radius) and the capture cross section 
for ions and electrons. Therefore, it is important to adopt a 
realistic model for the porosity of aggregates. 

In th is study, we adopt the porosity model developed by 
IOTS091 This model is based on iV-body simulations of suc- 
cessive collisions between aggregates of various sizes. This 
model provides a natural ext ension o f the classical hit-and- 
stick aggregation models (see OTS09 and references therein). 
Collisional fragmentation and restructuring is not taken into 
account, so the porosity increase only depends on the phys- 
ical sizes of colliding aggregates. This assumption is valid 
as long as the collisional energy is sufficiently lower than the 
critical energy for the onset of collisional compaction. The 
validity of this assumption will be discussed in Section 5.4. 

2.3.1. Porosity Increase After Collision 
Our porosity model measures the size of a porous aggregate 
with the characteristic radius a = [(5 /3N)J2H =l (x.k-X) 2 ] 1 / 2 , 
where N is the number of constituent monomers within the ag- 
gregate, is the coordinate of the k-th constituent monomer, 
and X is the center of mass. Figure[TJshows the characteristic 
radius as well as the projected area A of a numerically created 
porous aggregate. In our model, the porosity of each aggre- 
gate is characterized by a and N, while the projected area A 
is assumed to be a function of them. In the following sub- 
sections, we summarize how a and A are calculated in this 
model. 

The porosity evolution of aggregates after a collision is ex- 
pressed in terms of the increase in the porous volume V = 
(47r/3)a 3 . For a collision between aggregates with volumes 
V\ and V 2 (^ Vi), the volume of the resulting aggregate, V1+2, 
can be generally written as 

Vi + 2 = Vi + (l + x)V 2 , (17) 

where x is a dimensionless factor depending on V\ and V2. We 
refer to \ as the "void factor" since it identically vanishes for 
compact aggregation. 

It is known that ther e are two limiting cases forhit-and-stic k 
collisions (see, e.g.. iMukai et alj|1992t iKozasa et al.lll993l) . 
One is called the ballistic cluster-cluster aggregation (BCCA) 
where aggregates grow only through equal-sized collisions. 



On average, the characteristic radius of BCCA clusters is re- 
lated to the monomer number N as 



obcca ~ o.qN 



l/^BCCA 



(18) 



where a,Q is the radius of monomers and Dbc ca ~ 1-9 is 
the fr actal dimension of BCCA clusters (e.g., IMukai et al.l 
1992). The void factor for the BCCA growth can be calculated 
from Equation (QJi as X bcca = 2 3 /° BCCA -2 « 0.99 (lOTS09h . 
The opposite limit is called the ballistic particle-cluster ag- 
gregation (BPCA), in which an aggregate grows by collid- 
ing with monomers. On average, the characteristic radius of 
BPCA clusters is given by «bpca ~ (1 - PbpcaT 1 ^ 3 a^N , where 
Pbpca = 1 - ( AWabpca) 3 ~ 0-87 4 is the porosity of BPCA 
clusters (e.g.. IKozasa et al]|1993l) . The void fact or is found 
to be xbpca = Pbka/U -Pbka) « 6.94 (lOTS09h . Note that 
both xbcca and Xbpca are constant. 

To obtain x for more general cases, OTS09 presented a new 
aggregation model called the "quasi-BCCA" (QBCCA). In 
the QBCCA, an aggregate grows through unequal-sized colli- 
sions with a fixed mass ratio N2/N1, where Ni and A^(< N\) 
are the monomer numbers of the target and projectile, respec- 
tively. The projectile is chosen among the outcomes of earlier 
collisions, so that the resultant aggregate has a self-similar 
structure. OTS09 performed A^-body simulations of aggre- 
gate collisions with various size ratios and found that the void 
factor for QBCCA is approximately given by 

Xqbcca(Vi/V 2 ) = Xbcca- 1.031n(— ^— y) . (19) 

Note that xqbcca approaches to xbcca in the BCCA limit 
(V1/V2 -> 1) as must be by the definition of BCCA. 

Unfortunately, Equation f!9[ does not reproduce the void 
factor in the BPCA limit (Vi /V2 —> 00) . To b ridge the gap be- 
tween the BCCA and BPCA limit, IOTS09I considered a for- 
mula 

X = min{xQBCCA(Vi/y 2 ), Xbpca}- (20) 

It is easy to check that Eqaution ( f20b approaches to xbcca and 
Xbpca in the BCCA and BPCA limits, respectively. Equa- 
tion d20b will be used in the numerical simulations presented 
in Section 4 to determine the porosity (volume) of aggregates 
after collisions. 

2.3.2. Projected Area 

The projected area A is another key property of porous ag- 
gregates. This does not affect only the charge state of the 
gas-dust mixture (Equation (0]i) but also the drift velocity of 
individual aggregates (Equation ([Toll). 

For BCCA clusters, th e projected area ave raged for fixed N 
is well approximated by (Mina toet alJ l2006) 



AfiCCA = TTfln x 



12.5Af ' 685 exp(-2.53/A^ 0920 ), N< 16, 

0.352A^+0.566Af a862 , N ^ 16. 

(21) 

For BPCA clusters, the averaged projected area is simply 
related to the radius as A B pca ~ ™ 2 - For more general 
porous aggregates, including QBCCA clusters, the averaged 
projecte d area is well approximated by (Equation (47) of 
(OTS03)) 



A = 



1 



1 



1 



^bcca(A0 710 ira BC CA(N) 2 



(22) 
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where a and N are is the characteristic radius and monomer 
number of the aggregate considered, and a^ccpSN) and 
Abcca(^V) are the characteristic radius and projected area of 
BCCA clusters with the same monomer number Af (i.e., Equa- 
tions ( f]~8T > and d2"Tl) ). respectively. Note that the above formula 
reduces to Equation ( l2Tb in the BCCA limit (a s» abcca) and 
to A s» ira 2 in the BPCA limit (a <C «bcca, ™ 2 "C Abcca)- 

It should be noted that the above formulae can be only used 
for the average value of A. This does not bother us when we 
compute the charge state of aggregates, since it only depends 
on the total projected area A tot . However, we cannot ignore 
the dispersion of A when we calculate the differential drift 
velocity between aggregates, especially between BCCA-like 
clusters. For example, let us consider two BCCA clusters with 
different masses Ni and A^^iVi). As Equation (ETT i suggests, 
the mean mass-to-area ratio N /Abcca of BCCA clusters ap- 
proaches to a constant value in the limit of large N. Hence, if 
we ignored the area dispersion, we would have a differential 
drift velocity Auo oc A{N/A) vanishing for very large N\ and 
N 2 even ifN\ • Clearly, this would lead to underestimation 
of Aud and overestimation of the electrostatic repulsion. 

To avoid this problem, we should replace \N\ /A\ -N 2 /A 2 \ 2 

in £ D with \Ni/A l -N 2 /A 2 \ 2 , not with \N\/M - N 2 /T 2 \ 2 , 
where the overlines denote the statistical average. In partic- 
ular, if the standard deviation of N/A scales linearly with its 
mean, we can write [A(N/A)] 2 as (see Appendix) 



[A(N/A)]=\N l /A l -N 2 /A 2 \ +e 2 Y,(Nj/Ajy, (23) 

7=1,2 

where e is the ratio of the standard deviation to the mean of 
N/A. In the Appendix, we evaluate e from the numerical data 
on the projected area of sample BCCA clusters. We find that 
e can be well approximated as ~ 0.1 for N < 10 . In the fol- 
lowing sections, we will assume e = 0.1 for all aggregates, 
since the area dispersion is only important for collision be- 
tween BCCA-like clusters. 

2.4. Nondimensionalization 

As seen above, our dust model is characterized by a number 
of model parameters. To find a truly independent set of model 
parameters, we scale all the physical quantities involved into 
dimensionless ones. 

We introduce the dimensionless radius and mean projected 
area, 



7Tflg 



(24) 



(25) 



Also, we scale the mass M with the the monomer number N = 
M/mo, where mo is the mass of monomers. The normalized 
drift energy £ D and electrostatic energy £ E are already given 
by Equations (TBI) and ( fTBT ). respectively. Using {1Z, A, N) 
instead of (a, A, M), we have 



£d - Sd 



N1N2 
N1+N2 



N2 
A 2 



+ 6 



2 E 

7=1,2 



n x n 2 

TZi+TZ 2 : 



(26) 



(27) 



where the dimensionless coefficients fo and fg are defined as 

2 



Si 



D - 



m / gp a Q / nm g 



2k B T V Pg V V 



1.7 x 10 



-5 ( fl o A 5 1 Po x 3 



0.1 fj,mJ V 1 g cm 3 



)(: 



. 10 3 cm s 2 / \ 10 11 s cm 



-2/ T \" 2 

(took) (28) 



with the monomer material density po = 3mo/47TflQ. 
We also introduce the normalized distribution function 

n(M)dM 



F(N)dN : 



(30) 



where «o is the number density of monomers in the initial 
state. Note that the mass conservation ensures J NJ-(N)dN = 
1 . Using T, we rewrite the ionization parameter O as 



9: 



-4-totCtot 



(31) 



where A M = J A(N)T(N)dN and C tot = / K(N)T(N)dN are 
the normalized total projected area and capacitance, and h is 
a dimensionless ionization rate defined by 



h: 



C,n g e 2 



8.1 x lO" 3 *^, 



Pg 



ao 



Po 



( Pd/Pg Y 
V 0.01 J 



do- 11 



0.1 iimJ V 1 g cm 3 
1 / T \ _3 / 2 / C 



gem 



100 



k)" (10^) " (32) 



The surface potential is determined as a function of 8 by 
Equation (O, or 

1 expCE-Wpo) | AotCtot g =Q „ 



1 + * 



l + * c 



h * r 



where we have eliminated SjUj/s e u e using Equation ©. 

From the above scaling, we find the collisional growth of 
charged dust aggregates can be characterized by five dimen- 
sionless parameters (fo, /e, h, e, ^oo). 

3. MONODISPERSE GROWTH MODEL 

Before proceeding to the full simulations, we consider sim- 
plified situations where dust grows into monodisperse aggre- 
gates, i.e., where all the aggregates have the same monomer 
number Af at each moment. This greatly helps us to under- 
stand the results of the numerical simulations shown in the 
following section. 

Within the framework of the hit-and-stick aggregation 
model, the monodisperse growth is equivalent to the BCCA 
growth. Thus, the assumption of the monodisperse growth is 
expressed by the following relations: 



a = a ( — 

V mo / 



K = N" D , 



A = Abcca (AO 



A = A(N) : 



l BCCA 



(AO 



(34) 



(35) 
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„(M')= — 5(M'- 
M 



At) ^ T(N') = ^S(N-N'), (36) 



where D is the fractal dimension of BCCA clusters and S(x) 
is the delta function. Since D is close to 2 (see Section 2.3.1), 
we simply set D = 2 in the following calculation. Note that the 
l/N factor appearing in Equation (l36*l l accounts for the mass 
conservation J NT(N)dN = 1. 

Under the monodisperse approximation, the drift and elec- 
trostatic energies (£o and £e) can be given as a function of N. 
Substituting Equations d34t and (l35l l into Equation (|26| |. the 
drift energy can be written as 



So = foe' 



/V 3 



A(N) 2 



(37) 



Thus, under the monodisperse approximation, /b and e de- 
generate into a single parameter /be 2 . Similarly, the electro- 
static energy is written as £e = (/£ , /2)(^'/^ r 00 ) 2 A^ 1 / 2 , where 
* is given by Equation ([33]l with A tot = A(N)/N and C m = 
1Z/N = N~ x l 2 . The expression for £e can be further simpli- 
fied using the approximate formula for ^ (Equation (JVj) to 
eliminate *f> /^oq. The result is 



c fE 

OF = -=- 



1+ h 



N V2 



-0.8 



A(N)J 



-2.5 



N 1 ' 2 . 



(38) 



Note that this expression no longer involves ^oq. From Equa- 
tions d3~7T > and d38l . we find that the outcome of the monodis- 
perse growth is (approximately) determined by three parame- 
ters /be 2 , /#, and h. 

For later convenience, we define the "effective kinetic en- 
ergy" £ K as 

£k=1+£d, (39) 

or equivalently, Ek = £k^T = k%T + M ^(Auo) 2 /2. The first 
term in the right hand side of Equation (l39l accounts for the 
contribution of Brownian motion to the collisional energy (~ 
k^T). We expect that the monodisperse growth is strongly 
suppressed when £e exceeds £k- 

Here, we give some examples to show how £k and £e de- 
pends on the parameters. Figure|4]shows £k as a function of N 
for foe 2 = 10~ 7 . As found from this figure, the kinetic energy 
is constant alN< 10 6 due to Brownian motion (£% and 
increases with mass at N > 10 6 due to the differential drift 
(£k ~ £d N 3 /A 2 ) . The qualitative behavior is the same 
for every /oe 2 . The value of foe 2 only determines the mass 
at which the differential drift starts to dominate over Brown- 
ian motion in the kinetic energy. In figure @] we also plot £e 
for/ £ = 10 with varying the value of h(= 10" 4 5 , 1(T 6 , 1(T 7 - 5 ). 
For all the cases, £e quickly increases with N and finally be- 
comes proportional toTZ = N^ 2 . This reflects the transition of 
the plasma state from the IDP (* rj 9 oc N 3/2 /A) to the IEP 
(<]/ w ^oo). In the IEP limit, £e depends on /e but is inde- 
pendent of h. An important difference among the three exam- 
ples is the timing of the plasma transition: for smaller h, £e 
approaches the IEP limit at larger A^. This difference makes 
the ratio between £e and £ # qualitatively different among the 
three cases. For h = 10~ 4 5 , £e exceeds £k when the relative 
motion is dominated by Brownian motion. For h = 10~ 6 , by 
contrast, £e exceeds £k when the relative motion is dominated 
by the differential drift. For h = 10~ 7 5 , £e does not exceed £k 
for arbitrary A^. As we see in Section 4, this difference is a 




10 3 10 6 10 9 10 12 
mass N = M/m 

FIG. 4. — Examples of the effective kinetic energy Ek = 1 +£d and the 
electrostatic energy £e as a function of N. The black thick curve shows Ek 
for /fit 2 = 10~ 7 , and the three gray curves show £e for /e = 10 and h = 10~ 4,5 , 
1CT 6 , and 1(T 7 - 5 . The black arrow shows the critical drift mass Nd defined 
in Section 3.1, while the gray crosses show the freezeout mass Nf defined in 
Section 3.4 for h = 10^ 5 and lO" 6 . For h = 1CT 7 5 , £ E is below £ K for all N, 
so the freezeout mass is not defined. 

key to understand the collisional growth of dust aggregates 
with size distribution. 

To quantify these differences for general cases, we intro- 
duce the following quantities: 

• The drift mass No- This is defined as the mass at which 
the relative motion starts to be dominated by the differ- 
ential drift. 

• The plasma transition mass Np. This is defined as the 
mass at which the plasma state shifts from the IDP to 
the IEP. 

• The maximum energy ratio (£e / £ic)max- This is the 
maximum value of the ratio £e/£k in the monodisperse 
growth. If (^/f(;) max > 1, the electrostatic energy £e 
exceeds the kinetic energy £k at a certain mass. 

• The freezeout mass Nf- This is the mass at which £e 
starts to exceed £k- Note that the freezeout mass is only 
defined when (£ E / £ A:) m ax > 1 ■ 

In the following subsection, we describe how these quantities 
are related to the parameters (/be 2 , /£, h). 

3.1. N D : the Drift Mass 

The first and second terms in the right hand side of Equa- 
tion ( l39l represents Brownian motion and the differential 
drift. Since the second term monotonically increases with 
Af, there exists a critical mass at which the dominant relative 
motion changes from the Brownian motion to the differential 
drift. We define Nd as the critical mass satisfying £e,{Nd) = 1. 
Using Equation ( 137) . the equation for No is written as 



A(N D ) 2 



■foe 2 



(40) 



This equation implicitly determines No as a function of /be 2 . 
For example, No ~ 3 x 10 6 when foe 2 = 10~ 7 (see Figure|4]l. 

Figure[5]shows the solution to Equation ( l40b - When foe 2 <C 
1 , Nd is well approximated as 



A'd : 



1 



b 2 f D e 2 



(41) 
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contours: N D (solid), N P (dashed) 
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FIG. 5. — Contour plot of the drift mass Np (Equation )40t ; solid lines) and 
the plasma transition mass Np (Equation 143 1 ; dashed lines) as a function of 
h (x-axis) and foe~ (y-axis). 




log N 



log N 



FIG. 6. — Schematic diagrams describing the mass dependence of the effec- 
tive kinetic energy Eg (a) and the electrostatic energy Ee (b). Here, No and 
Np are t he drift mass and plasma transition mass defined by Equations i40\ 
and )43t . respectively. The dashed lines with arrows indicate how Ex and Ee 
depends on the parameters /be 2 , fs, and h. 



where b = 1 /0.352 = 2.84 is the mass-to-area ratio N/A(N) 
in the limit of N —> oo. Using Equation fiTT i. £k is simply 
rewritten as 

N 

£k~1 + —, (42) 

which asymptotically behaves as £k ~ 1 for N -C No and 
£ K w N/Nd for AT ^> Afo. The asymptotic form of is 
schematically illustrated in Figure [6f a). 

3.2. A^>: /fte Plasma Transition Mass 

Another important quantity is the critical mass at which the 
plasma state changes from IDP to IEP. We define the critical 
mass Np such that <d(Np) = (see Equation ©). Using 




h 

FIG. 7. — Contour plot of the maximum energy ratio (Ee/Ek^hx divided 
by /e as a function of h (x-axis) and foe 2 (y-axis). The dashed line represents 
No = Np (see also Figure [5j- The two parameter regions (I) and (II) are 
characterized by No ^S> Np and No -C Np, respectively (see also Figure[8). 



Equation QTl ). this condition can be written as 
A(N P ) 



n: 



3/2 



■ h. 



(43) 



Note that Np depends on h only. 

Figure [5] shows the solution to Equation d43l as a function 
of h. If h <C 1, Np is well approximated as 



N P : 



b 2 h 2 ' 



In this case, £e can be approximately written as 

-2.5 



C ~ f E 
OF 



1 + 



1/2 



(44) 



(45) 



which asymptotically behaves as £ E ~ (fE/2)N 3 / 2 /Np for 
N < Afr and as £ £ w (f E /2)N 1 ' 2 for AT > Afc. The asymp- 
totic form of £ £ is illustrated in Figure[6fb). 

3.3. (£ e/£ A:)max-' the Maximum Energy Ratio 

The maximum energy ratio (£ £ / £A:) ma x determines whether 
the electrostatic energy exceeds the kinetic energy during the 
monodisperse growth. Since £e scales linearly with / £ , the 
quantity 1 (£ £ / £ K)max depends only on foe 2 and h. 

Figure [7] plots /e 1 (£e /£i()max as a function of foe 2 and 
ft. It is seen that the maximum energy ratio behaves differ- 
ently across the line No = Np. This can be easily understood 
from Figure [8] which schematically illustrates the mass de- 
pendence of £k and £e (Equations d42l and (05])). If No ^> 
Np (/dc 2 <C h 2 ), the energy ratio £ e /£k reaches the maxi- 
mum at N « AT D . Since £d(Nd) = 2 and £ £ (Afo) « f E N D ' 2 /2, 
we obtain (£ E /£ K )m^ ~ f E N D /2 /4 f E /Abf x D ,2 e indepen- 
dently of ft. If Ad < AV (/ D e 2 > ft 2 ), by contrast, £ £ /£/t 
reaches the maximum at AT s» A^>. Using £o(Np) = Np/No 
and £ £ (AW « / £ /2 3 5 , we have (£ E /£ K )m^ « f E N D /2 35 N P « 
f E h/2 35 bf D e 2 , which depends on both / D e 2 and ft. 
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log £ 




FIG. 8. — Schematic diagrams describing the dependence of the maximum 
energy ratio (£e I £K)mnx on No and Np in regions (I) and (II) shown in Fig- 
ure[7| The black and gra y lines shows the asymptotic behavior of Ek and Ee 
(Equations (42) and 145) ) as a function of N, respectively. When No 3> Np, 
or equivalently fof. <C h (region (I); upper panel), the energy ratio max- 
imizes at N ~ No- In the opposite limit (region (II); lower panel), £e/£k 
maximizes atNzzNp. 



contours: N F (/e=10) 




1CT 8 10" 6 10" 4 1(T 2 10° 10 2 
h 



FIG. 9. — Contour plot of the freezeout mass Nf (thin solid curves) for 
/e = 10 as a function of h (x-axis) and foe 2 (y-axis). The dashed and dotted 
curves show £e(Nd) = 2 and £e(Np) = 1, respectively. The regions (i), (ii), 
an d (hi ) are characterized by the values of £e(Nd) and £e(Np) (see also Fig- 
ure llOt . Above the thick solid curve (region (iv)), the maximum energy ratio 
(Ee / EKjaaa. is less than unity, so the freezeout mass is not defined. 
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FIG. 10. — Schematic diagrams describing the location of the freezeout 
mass Nf in the mass space for three parameter regions (i), (ii) and (iii) 
shown in Figure [9] The blac k an d gr ay li nes shows the asymptotic behav- 
ior of Ek and £e (Equations I42t and 145 1 ) as a function of N, respectively. 
If £e(Nd) 2> 1 (cases (i) and (ii); top and middle panels), £e exceeds Ek 
in the Brownian motion regime (i.e., Nf <S No)- If £e(Nd) <C 1 but still 
(£e / £x)max 3> 1 (case (iii); bottom panel), £e exceeds £k in the differential 
drift regime (i.e., Nf 3> No). 



/e = 10. We see that Nf depends on these parameters dif- 
ferently depending on the values of £e(Np) and £e(Nd). To 
understand this, in Figure [10] we schematically show Ek and 
Ee as a function of for the three cases. If £e(Nd) ^> 1, £e 
starts to exceed £k when the relative velocity is dominated by 
Brownian motion (i.e., Np -C No)- In this case, the condition 
determining Nf is given by £e(Nf) ~ 1, which implies Nf ~ 
(2// £ ) 2 for £ (Np) « 1 and N F « {2N P /f E ) 2 ^ ps (2/b 2 f E h 2 ) 2 / 3 
for £(N P ) > 1. If £ E (N D ) < 1 but still {£e/£k)^ > 1, £ E 
exceeds £k after the relative velocity is dominated by the dif- 
ferential drift (i.e., Np 3> No). In this case, the condition for 

N F is given by ( f E /2)N^ 2 /N P s» N F /N D , hence Np is given 
by Np w (2N P /f E N D ) 2 w (2f D e 2 /f E h 2 ) 2 . 



3.4. iVf/ the Freezeout Mass 

When (Ee/Ek)™™ > 1, there exists a critical mass Afc at 
which the electrostatic energy £e takes over the kinetic energy 
£ K . As we will see in Section 3.5, the monodisperse growth 
is strongly suppressed at > Nf. For this reason, we refer 
to Nf as the "freezeout mass." The freezeout mass can be 
calculated from the condition £k(Nf) = £e(Nf) once the three 
parameters /^e 2 , fp, and h are specified. 

In Figure [9] we plot Nf as a function of f D e 2 and h for 



3.5. The Outcomes of Monodisperse Growth 

As mentioned above, the monodisperse growth is ex- 
pected to slow down at the freezeout mass N w Nf when 
(£e/ '&)max > 1. Here, we demonstrate this by numerically 
calculating the mass evolution. 

Under the monodisperse approximation, the evolution of 
aggregate mass is given by 



dM 



dN 

df =>C 



(46) 
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TABLE 1 

Critical Masses and the Maximum 
Energy Ratio for (Jd^Je) = (10^ 7 , 10) 



h 


N D 


N P 


(£e 1 '^A:)max 


N F 




10 6.3 


10 7 - 2 


10 20 


10 51 


icr 6 


10 6.3 


10 io.i 


10 0.5 


10 88 


io- 7 - 5 


10 6.3 


10 13.1 


10 -i.o 






time T 

FIG. 11. — Mass evolution in the monodisperse model calculated from 
Equation |46j for f D e 2 = lfr 7 and f E = 10 with various values of h. The black 
arrow indicate the drift mass No. while the lower and upper arrows show the 
freezeout mass Nf for h = 10 -4 ' 5 and 10~ 6 , respectively. The evolution for 
the uncharged case (i.e., h = 0) is shown by the dashed curve. 

where T = nQira^t J 8£ B T /ttihq and JC = 

K I (ttciq \J%k-&T I TTino) are the dimensionless time and colli - 
sional rate coefficient. We numerically solve Equation d46l > 
with initial condition N(T = 0) = 1 . 

As in the beginning of this section, we consider three cases 
of h = 10" 4 5 , 10" 6 , and 10" 7 5 with fixed f D e 2 = 10" 7 and f E = 
10. Listed in Table[T|are the critical masses (No, Np, Nf) and 
the maximum energy ratio {£ E /£K)ms.\ for these cases. We 
also consider the uncharged case with the same value of /be 2 . 

3.5.1. Without Charging 

In Figure QT| the mass evolution for the uncharged case is 
shown by the dashed curve. The black arrow in the figure in- 
dicates the critical drift mass No = 10 6,3 . We find that the mass 
grows as T 2 until reaching No, and then grows exponentially 
with T. This evolutionary trend can be directly proven from 
Equation (1461 . Without charging, the collision kernel JC is just 
the product of the geometrical cross section cx 1Z 2 = N and the 
relative velocity Am. When /V <C No, the relative velocity is 
dominated by Brownian motion (i.e., Am cx TV -1 / 2 ), and we 
have JC cx 1Z 2 N~ l l 2 cx TV 1 / 2 . Inserting this into Equation (146) . 
we have /V cx T 2 . When /V No, by contrast, the relative 
velocity is dominated by the differential drift (Am cx N / A), 
and hence JC cx N7Z 2 /A. Since the projected area A roughly 
scales with 1Z 2 , we have JC cx N. Hence, from Equation (l46t . 
we find N cx exp(fiT), where ft is a constant growth rate. 

3.5.2. With Charging 

The mass evolution for the charged cases is plotted in Fig- 
ureQT]by g ra y curves. The gray arrows in the figure indicate 



the freezeout mass Nf for h = 10 4,5 and 10 6 . As expected, 
we observe significant slowdown in the growth at N sa Nf for 
the two cases. At T= 10 4 , the aggregate mass is N ~ 10 5,7 for 
h = 10" 4 5 and N « 10 s 9 for h = 10~ 6 , which is consistent with 
the predicted freezeout mass (see Table [TJ. We have com- 
puted the mass evolution for the two cases until T = 10 6 , but 
the final masses 10 5 9 and 10 are not very different from the 
values at T = 10 4 . 

For h = 10~ 7 5 , by contrast, the evolution curve of N is in- 
distinguishable from that for the uncharged case, meaning that 
the electrostatic repulsion hardly affects the aggregate growth. 

To summarize, we have confirmed that dust can continue 
the monodisperse growth only if 

(t-) < 1- (47) 

V Ofc ' max 

4. NUMERICAL SIMULATIONS INCLUDING SIZE 
DISTRIBUTION 

As shown in the previous section, dust aggregates could not 
grow beyond the freezeout mass Nf if the condition Wh is not 
satisfied and if the size distribution were limited to monodis- 
perse ones. In this section, we study how the outcome of dust 
growth changes when we allow the size distribution to freely 
evolve. 

To compute the evolution of size distribution, we em- 
ploy th e "extended" Smoluchowski method developed in 
OTS09. This method treats the number density n(M) 
and the mean volume V(M) of aggregates with mass M 
as time-dependent quantities, and calculates their tempo- 
ral evolution simultaneously. This method allows us to 
follow the porosity evolution consistently with collisional 
growth, which cannot b e done with the conven t ional Smolu- 
chowski method (e.g.. iNakagawa et alj I198U iTanaka et alj 
l200HlDullemond & Dominikll2005h . 

In the extended Smoluchowski method, the temporal evo- 
lution of n(M) and V(M) is given by two equations, 

dn(M) _ 1 r dM >K( M r. M _ M >) n ( M >) n ( M _ M ^ 
dt 2 J 

-n(M) I dM' ~K(M;M')n(M'), (48) 
Jo 

d[V(M)n(M)] 1 f M , - . .- . 

' = - / dM' V 1+2 (M';M-M')K(M';M-M') 
dt 2 J 

xn(M')n(M-M') 

-V(M)n{M) / dM'K(M;M')n(M'), (49) 
Jo 

where K(M\;M2) and Vi+2(Mi;M2) are the collisional rate co- 
efficient K (Equation ( fT2l ) and the aggregate volume V 1+ 2 af- 
ter a collision (Equation ( fTTI i) evaluated for V\ = V{M\) and 
V2 = V(M2). In this study, we determine V\+2 using the for- 
mula for hit-and-stick collisions (Equation d20li). 

We numerically solve Equations d48l ) and d49b using the 
fixed bin scheme described in 1OTS091 This scheme divides 
the low-mass region mo ^ M ^ Afbdino into linearly spaced 
bins with representative masses = kmo (k= 1,2, .. . ,J\fbd) 
and the high-mass region M > Afbd^o into logarithmically 
spaced bins with M k = 10 1 l Mu M%-\ (k = J\f bd + 1, . . .). The 
number Afbd controls the resolution in the mass coordinate. In 
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FIG. 12. — Evolution of the mass distribution function J-'iN) (upper panel) 
and the mass-radius relation TZ(N) (lower panel) for the uncharged case of 
(/be 2 ,e) = (10~ 7 , 10 _1 ). The gray curves show the snapshots ofiV 2 J r (iV) and 
Tl(N) at various times, T = 10 1 , 10 1 5 , 10 2 , . . . , 10 4 (from left to right). Note 
that the curves for H(N) overlap each other. The arrows indicat e th e criti- 
cal mass No calculated from the monodisperse theory (Equation j40t ). The 
crosses an d op en circles in the upper panel indicate the averaged mas s (N) 
(Equation )50» and the weighted averaged mass {N) m (Equation j5U ), re- 
spectively. In the lower panel, the mass-radius relations for the fractal dimen- 
sions of D = 2 and 3 are shown by the dashed and dotted curves, respectively. 



this study, we set Nbd = 80 (meaning Mk+\ /M^ = 1 .03 for the 
high-mass range). The temporal evolution is computed using 
the explicit, forth-order Runge-Kutta method. The time in- 
crement Af for each time step is continuously adjusted so that 
the fractional decrease in the number density during Af does 
not exceed 8 t for all bins, where 8, is a constant parameter. 
We take 8, = 0.02 in the following calculations. 

4.1. Without Charging 

Figure [T2] shows the solution to Equations d48l and d49l 
for the uncharged case of (/oe 2 , e) = (10~ 7 , 10 _1 ). The upper 
panel displays the mass distribution function J-(N) at various 
times T. Note that the vertical axis of this panel is chosen 
to be N 2 F{N), which is proportional to the mass density of 
aggregates belonging to each logarithmic mass bin. 

To characterize the evolution of the mass distribution, we 
introduce the average mass (Af) and the mass-weighted aver- 
age mass (N) m defined by 



1 



J™NT(N)dN 



N \F(N)dN, (51) 



J™F(N)dN / J"(A0rfJV' 



(50) 



where we have used the mass conservation L NJ-(N)dN = 1 . 
Note that (Af) is inversely proportional to the total number 
density of aggregates, /„ T(N)dN. Roughly speaking, (N) 
represents the mass scale dominating the number of aggre- 
gates in the system, while (N) m represents the mass scale 
dominating the mass of the system. Also note that (N) m 
can be written as (N 2 )/(N), and the dispersion (8N 2 ) = 
(N 2 ) - (N) 2 of the mass distribution is written as (8N 2 ) = 
(N) 2 ((N) m /(N) - 1). Hence, the ratio (N) m /(N) measures 
how the mass distribution deviates from the monodisperse dis- 
tribution. In the upper panel of Figure [T2] we indicate (N) 
and (N) m at each time with crosses (x) and circles (o), re- 
spectively. 

The evolution of the mass distribution can be divided into 
two stages. During (A^),,, < No, the mass distribution evolves 
with small dispersion ((N) sa (N),„). The average masses (N) 
and (N)„, grow approximately as T 2 , which is consistent with 
the prediction of the monodisperse theory (see Section 3.5.1). 
These imply that the monodisperse approximation is good 
when Brownian motion dominates the relative motion of ag- 
gregates. 

However, the monodisperse approximation becomes less 
good once (N) m exceeds No- In this stage, we observe a 
power-law tail extending from N sa (N) m down to N « A^d. In 
fact, we see that the growth rate of (N) m (i.e., dln(N) m /dT) 
is approximately twice as high as that of (N). This means 
that the relative width of the distribution ((SN 2 ) 1/2 /{N) = 
y/ (N)„,/ (N) - 1 « y/ (N) m / (N)) increases exponentially with 
timfl As we will see in the following subsection, the broad- 
ening of the mass distribution plays a key role when dust 
charging is present. 

The lower panel of Figure Q~2] shows the temporal evolu- 
tion of the mass-radius relation 7Z(N). We see that 7Z(N) 
approximately obeys a fractal relation TZ « N l ' D , where the 
fractal dimension is D w 2 independently of the time (see 
the dashed line in the panel which shows the exact rela- 
tion TZ = N x/2 ). This fact validates the assumption TZ = N l/2 
made in the monodisperse theory (see Section 3). In fact, 
the fractal dimension close to 2 is a general consequence of 
dust growth without collisional compaction when aggregate 
collision is driven by Brownian motion and differential drift 
(OTS09). Detailed inspection shows that values D = 1.95 and 
2.03 better fit to the data if the fitted region is limited to the 
Brownian motion regime (N < No) and the differential drift 
regime {N > No), respectively. The differential drift leads to 
a slightly higher fractal dimension than Brownian motion be- 
cause the former reduces th e collis ion rate for similar-sized 
aggregates (see Figure 15 of lOTS09h . 

4.2. With Charging 

Now we show how the charging alters the evolution of the 
size distribution. As in Section 3, we consider three cases 
of h = 10" 45 , 10" 6 , and 10" 7 5 with (f D e 2 , f E , e) = (10" 7 , 10, 
10" 1 ). 

7 As pointed out by the referee, this is a general consequence of the kernel 
K. scaling linearly with the masses of colliding aggregates (this is the case for 
our kernel atW> No, see Section 3.5.1). In fact, the growth rate of (N),„ 
is known to be exactly twice as high as that of (Af) when the kernel is o f the 
form K(Ni ;N 2 ) oc N\ +N 2 (see, e.g., Figure 1 of[prmel & Spaans 20081). 
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FIG. 13. — Same as the upper panel of Figure [TSl but for three charged 
cases, h = 1CT 4 ' 5 , 1(T 6 , and 10~ 7 5 (from top to bottom). The other parameters 
are set to (f D e 2 ,f E , e, *) = (KT 7 , 10, 10 _1 , 10° 5 ). The gray arrows indicate 
the freezeout mass Nf predicted from the monodisperse theory. The dotted 
curves in the middle and bottom panels show the mass distr ibut ion when the 
surface potential ^ exceeds the critical value (Equation fl53» . 



In Figure [13] we show the temporal evolution of the mass 
distribution J-(N) for the three cases. The mass-radius re- 
lation 1Z(N) is not shown here because it is very similar to 
that for the uncharged case. For h = 10~ 45 , the monodisperse 
theory gives {£e / Zidmax > 1, predicting the freezeout of the 
growth at N « Nf ps 10 5 2 (see TableHJ. As expected, the evo- 
lution of the mass distribution starts to slow down at N « Nf, 
ending up with nearly monodisperse distribution peaked at 
N ~ 10 6 . In the simulation, we have followed the evolution at 
7~= 10 6 , but observed no significant growth after T > 10 4 . 

For h = 10~ 6 and 10~ 7,5 , by contrast, the outcome is qualita- 
tively different from the prediction by the monodisperse the- 
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FIG. 14. — Evolution of the average mass (N) (upper panel) and the 
weighted average mass (N) m (lower panel) as a function of time T. The 
gray curves indicate the results for three charged cases of h= 1CT 4 ' 5 , lfr 6 , 
10~ 7 , while the black dashed curves are for the uncharged case (h = 0). The 
other parameters are set to (/d./s, e, *oo) = (1CT 5 , 10, 10~' , 10° 5 ). The gray 
and black arrows indicate the critical drift mass No and the freezeout mass 
Nf predicted by the monodisperse theory, respectively. 



ory, as is shown in the middle and bottom panels of Figure[T3l 
respectively. For the case of h = 10~ 6 , the prediction was that 
the freezeout occurs at « Np ~ 10 9 . However, the simula- 
tion shows the size distribution evolving into a bimodal dis- 
tribution, in which one peak stays at N ~ No and the other 
continues growing towards larger N. Interestingly, similar be- 
havior is seen in the case of h = 10 despite the fact that the 
charging did not affect dust growth for this case within the 
monodisperse theory. 

The evolution of the size dispersion can be better under- 
stood if we look at the evolution of (N) and (N),„. Figure [141 
compares them among the three charged cases together with 
the uncharged case. See also Figure [TT] in which the predic- 
tion from the monodisperse theory is shown. For h = 10~ 4 5 , 
both (N) and (N) m evolves as the monodisperse theory pre- 
dicts. However, for h = 10~ 6 and 10~ 7 5 , (N) stops growing at 
certain values, while (N) m continues growing as for the un- 
charged case. This means that, in the latter cases, only a small 
number of aggregates continue growing but nevertheless carry 
the greater part of dust mass in the system. 

As we explain below, the transition to the bimodal distribu- 
tion can be characterized by three steps: 

1 . At (N} m > No, a long tail is formed at the low-mass end 
of the size distribution. 
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2. Since aggregates belonging to the low-mass tail have 
a relatively small kinetic energy, they stop growing 
as the surface potential '5 reaches a certain value '5* 
(see Equation (l53l below). These "frozen" aggregates 
provide the total capacitance C tot which no longer de- 
creases with time. This leads to the surface potential ^ 
of all aggregates no longer increasing with time. 

3. Consequently, aggregates of higher mass are less 
charged than in the case of the monodisperse growth. 
The growth of the high-mass aggregates is no longer 
inhibited by the charge barrier. 

The first step was already discussed in the previous subsec- 
tion. Here, we explain how the second step follows after the 
development of the low-mass tail. Let us approximate the 
mass distribution at the end of the first stage into two sub- 
groups, one representing the high-mass side and the other 
representing the low-mass tail. We characterize them with 
masses N\ ^> Nd and N 2 ~ Nd- The number of the low-mass 
aggregates decreases through their mutual collisions ("2-2 
collision") and through sweep-up by the high-mass aggre- 
gates ("1-2 collision"). This leads to the decrease in the total 
capacitance C tot and, in turn, the increase in the surface poten- 
tial *. We now write the relative kinetic energies for 1-2 and 
2-2 collisions as £k,12 and £k,22- Using Equations d26| i, fl39l l, 
and (|40]) together with N/A « b and N\ > N 2 ~ Nd, these en- 
ergies are approximately evaluated as £k.\i ~ 1 +2A^/A^> « 3 
and £k.h ~ 1 +N2/ND ~ 2, respectively. Note that £k.\i is 
nearly independent of N\ because the reduced mass is deter- 
mined by smaller aggregates and because the drift velocity 
oc N\ / A\ is nearly constant at large N\ . Meanwhile, the elec- 
trostatic energies (Equation d27| i) for 1-2 and 2-2 collisions 

are written as £ EiU a f E (V /*oo) 2 ft 2 « M*/*oo) 2 N D /2 

and £ E>22 a / £ (*/* oc ) 2 ^ 2 /2 « / £ (*/* 0O ) 2 ^ /2 /2, respec- 
tively. Again, £^12 is independent of N\, because the reduced 
radius is determined by smaller aggregates. Thus, the energy 
ratios for 1-2 and 2-2 collisions are obtained as 



£e,u _ f E V 2 N D ' 
£k. 



1/2 



12 



3* 2 



£e,22 _ fE&No' 

£k. 



22 



4* 2 



(52) 



independently of N\. Both the energy ratios exceed unity 
when "J > where 



1/2 



**: 



J E N l '\ 



Mi 



!/2 



1/2 



(53) 
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The values of C tot ,* for the two cases are indicated in the lower 
right panel of Figure [151 

We are now able to explain why the high-mass aggregates 
can grow beyond N sa Np in the case of h = 10 -6 . First note 
that they can grow only through their mutual collisions ("1- 
1 collision") because 1-2 collisions have been already in- 
hibited. The relative kinetic energy and electrostatic energy 
for 1-1 collisions are now given by £km ~ 1 +N\/Nd and 



£e,u « C/k/2)(**/* c 
tion ( T53l , we obtain 



) 2 N\ 12 . Using N > Nd and Equa- 



Note that is independent of h. For f D e 2 = 10 -7 and f E = 10, 
we obtain sa 0.02*00. 

The above consideration suggests that the freezeout of the 
low-mass aggregates occurs when exceeds the critical value 
\I f *. To confirm this, in the upper panel of Figure[l5] we plot 
\I/ versus the average mass (N) for h = 10 -6 and 10 -7,5 . We 
see that the increase in (N) stops when exceeds 

It should be noted that the evolution of * is also slowed 
down for * > This is because the "frozen" small aggre- 
gates govern the total electric capacitance C tot of the system. 
Using * rj 8 « bh^ oc/Ctot (as is for the IDP limit), the total 
capacitance when * s» \]/* can be evaluated as 



fcM^o ^ (bf E ) l ' 2 h 

** ^2(/ D e 2 ) 1 /4' 



(54) 



Nn\ 1/2 



1/2 



\Ni 



< 1. 



(55) 



Thus, we find that the energy ratio decreases with mass, and 
therefore the growth of the high-mass aggregates is no longer 
inhibited by the charge barrier. This is essentially due to the 
frozen aggregates keeping the surface potential * nearly con- 
stant. Without the frozen aggregates, * would increase as 



1/2- 

iVj , and the electrostatic energy £ E) u <x N^ would take 
over £k,\ 1 oc N\ at a certain mass as in the monodisperse case. 
With the frozen aggregates, by contrast, £b,u increases only 

1 12 

as Ay , so cannot exceed £k.\\- It should be noted that the 
increasing kinetic energy will in reality cause collisional com- 
paction at some stage, but this effect is neglected in our cal- 
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culation. 

One might wonder why the freezeout of the entire mass dis- 
tribution occurs for h = 10 -4 5 . The key difference between 
the two cases h = 10~ 45 and h = 10~ 6 is the timing at which 
the electrostatic barrier becomes effective. In the former case, 
the charge barrier becomes effective when the relative motion 
between aggregates is dominated by Brownian motion (i.e., 
Nf < No). In this case, the aggregates cannot overcome the 
barrier even if VE' is kept constant, since the electrostatic en- 
ergy £ E oc ^ 2 N x ' 2 grows with mass while the kinetic energy 
£ k ~ 1 does not. In the latter case, by contrast, the charge 
barrier becomes effective after the relative motion has been 
already dominated by the differential drift (i.e., No < Nf). In 
this case, the kinetic energy £k oc N can surpass the electro- 
static energy if \& is kept constant. 

Finally, we remark that "J* can exceed "Joo when fo^lfW 
is sufficiently large (see Equation (l53ll). In reality, however, 
the surface potential does not grow larger than ^oo. For such 
cases, the energy ratios in Equation (l52l never exceed unity, 
so we expect that low-mass aggregates do not stop growing. 
We will confirm this expectation in the following subsection. 

4.3. The Growth Criteria 

The above examples suggest that the criterion (£e / £ Ar) m ax ^ 
1 for the monodisperse growth no longer applies when the 
evolution of the size distribution is taken into account. To ob- 
tain a working criterion, we have performed numerical simu- 
lations for various sets of parameters (/oe 2 ,/£,/i). 

Figure[T6]shows the parameter space considered in the sim- 
ulations. We have chosen various sets of parameters (/be 2 , 
/e, h) for which (£e / '£/f)max falls within the range 0.1 .. . 10 3 . 
We have set e = 10 -1 in all of the simulations. 

We find that the outcome of dust evolution can be classi- 
fied into three types in terms of the temporal evolution of (N) 
and (N) m . In the first type, we observe that both (N) and 
(N) m stop growing at N ~ Nf. The outcome is characterized 
by frozen aggregates with a nearly monodisperse distribution 
peaked at N s» Nf as seen in the top panel of Figure [13] as 
seen in the top panel of Figure Qj] We will refer to this type 
of growth outcome as the total freezeout. In the second type, 
we see that (N) stops growing at a certain value while (N) m 
continues growing. The outcome is a double-peaked size dis- 
tribution consisting of low-mass aggregates frozen at N « No 
and ever-growing high-mass aggregates, as seen in the middle 
and bottom panels of Figure [l~3] We will call this type the bi- 
modal growth. In the third type, we observe that both (N) and 
(N) m continue growing. The outcome is a single-peaked dis- 
tribution of ever-growing aggregates as is for uncharged cases 
(see Figure [T2l>. We will call this type the unimodal growth 
to emphasize that the size distribution is characterized by a 
single peak. 

The outcome of the growth for each set of parameters is dis- 
played in Figure [16] Here, the crosses (x), filled circles (•), 
and open circles (o) show the parameter sets for which we 
have observed the total freezeout, bimodal growth, and uni- 
modal growth, respectively. It is seen that the total freezeout 
occurs for small foe 2 and large h, while the unimodal growth 
occurs when foe 2 is small. 

First, we examine whether the total freezeout regime can 
be well represented by a criterion of the form (£e /£s:)max > 
constant as suggested by the monodisperse theory (see Equa- 
tion d47li). In Figure [16] we show a criterion (£f/£s:)max > 3 
with the solid curve. It is seen that this criterion applies well 



TABLE 2 

Three Outcomes of the Growth of 
Charged Dust 



Conditions 




Outcome 


£e(N d )>6 




Total freezeout 


£e(N d )<6 **<*oo 


/4 


Bimodal growth 


£e(N d )<6 **>*oo 


/4 


Unimodal growth 



at large f D e 2 while it overestimates the size of the freezeout 
region at smaller foe 2 . It is clear that such a type of criteria 
do not explain the condition for the total freezeout to occur. 

However, a criterion applicable for all parameter ranges can 
be obtained if we slightly modify Equation (07}. The point is 
that the total freezeout is observed only in the Brownian mo- 
tion regime, i.e., only when the freezeout mass Nf is smaller 
than the drift energy No- This fact suggests that the total 
freezeout does not occur if (£e / '£/t)max ^ 1 but £e(No) -C 1 
(this is the case for the parameter region (iii) in Figures [9] and 
[Tol l. This expectation motivates us to introduce another en- 
ergy ratio, 

£e(N d ) £ e (N d ) 

£M) = —' (56) 

where we have used the definition of No, i.e., £k(Nd) = 2. 
Note that £e(No)/2 is the maximum value of £e/£k in the 
Brownian motion regime because £e monotonically increases 
with N and £k ^ 2 at N ^ No- In Figure [16] we show the line 
£e(Nd) = 6 with the dashed curve. We see that the line rep- 
resents the boundary of the total freezeout regime very well. 
Thus, we conclude that the criterion for the total freezeout to 
occur is given by 

£e(N d )>6. (57) 

A simple criterion is also obtained for the boundary be- 
tween the bimodal and unimodal growth regimes. As men- 
tioned in Section 4.2, the bimodal growth occurs only if the 
critical surface potential (Equation d53l l) is lower than ^qq. 
In Figure [16] we show the line ^ = x I'oo/4 with the dotted 
curve, with the dashed curve. We find that the condition for 
the bimodal growth to occur instead of the unimodal growth 
is given by 

**<^f- (58) 

To summarize, the outcome of charged dust growth can be 
classified into three cases (Table [2|. If £e(Nd) > 6, all ag- 
gregates stops growing before the systematic drift dominates 
their relative velocities. The outcome is a nearly monodis- 
perse distribution of frozen aggregates with typical mass f=s 
Nf- If £e(Nd) < 6 and < a large number of aggre- 

gates stop growing, but the major part of dust mass within the 
system is carried by a small number of ever-growing aggre- 
gates. If £e(Nd) < 6 and ^ > x I'oo/4, all aggregates continue 
growing with a single-peaked size distribution. The second 
case includes situations where no aggregates could continue 
growing if the size distribution is limited to a monodisperse 
one. This means that size distribution must be taken into ac- 
count when we discuss how the charging of aggregates affects 
their collisional growth. 

5. DISCUSSION 
5.1. An Application to a Protoplanetary Disk Model 
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FIG. 16. — Outcome of numerical simulations for various parameters. The crosses (x) show the parameters for which both (N) and (N)„, stop growing at 
N ~ Nf ("total freezeout"). The filled circles (•) indicate the parameters for which (N) stops growing at N ~ Nd while (N),„ does not ("bimodal growth"). The 
open circies (o) indicate where both (N),„ and (N) continue growing with a single-peaked distribution ("unimodal growth"). The solid, dashed, and dotted lines 
show where (&/&)max = 3, Ee(Nd) = 6, and = <I/oo /4, respectively. 

The growth criteria derived in Section 4 are general in a 
sense that no protoplanetary disk model is specified. Al- 
though application to particular disk models is the subject of 
Paper II, we will show here one example of how to use the 
criteria. 

Here, w e adopt the mi nimum-mass solar nebular (MMSN) 
model of (Hayashi 1981). In this model, the gas temperature 
T and the Kepler rotational frequency fiic are given by T = 
280(r/l AU)-'/ 2 K and n K = (2tt/1 yr)(r/l AU)" 3 / 2 rad s" 1 , 
where r is the distance from the Sun. The gas density 
p g and the vertical component of the stellar gravity g are 
given by ^ = 1.4 x 10" 9 (r/l AU)- I1/4 exp(-z 2 /2// 2 ) g cm" 3 
and g = fllz = 0.020O/1 AV)- y4 (z/H), where z is the dis- 
tance from the midplane of the disk and H = c,/S1k = 5.0 x 
10 n O/l AU) 5 / 4 cm is the gas scale height. In this subsec- 
tion, we neglect the effect of disk turbulence to dust col- 
lision and assume the stellar gravity as the only source of 
dust differential drift. For the material density of monomers 
and the dust-to-gas mass ratio, we ignore the sublimation of 
ice for simplicity an d set po = 1.4 g cm" 3 and Pd/p g = 0.014 
(Tana kaet al.l l2005). The maximum surface potential is 
taken to be 2.81 as is for m, = 24wh and = 0.3. Substituting 
these relations into Equation d28l l. j29l , and (|32T > and setting 
z = H,we obtain 
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FIG. 17. — Map of the minimum-mass solar nebular model in the h-fo 
plane. The thick solid line shows the radial profile of h (x-axis) and fo Cy- 
ans) at one scale height above the midplane of the disk, with the filled squares 
indicating the distances from the central star. The break in the line approxi- 
mates attenuation of cosmic-rays and X-rays at inner radii. The gray region 
below the dashed curve indicates where we predict the to tal f reezeout of frac- 
tal dust growth (see the freezeout condition, Equation \51\ ). Note that we 
have used a relation between fg and fo to project the freezeout region onto 
in the h-fo plane (see text). The thin solid line shows (£e / £K )ma\ = 3; fractal 
dust growth beyond the electrostatic barrier is possible between this line and 
the dashed line because of the effect of dust size distribution (see Section 4). 



/, = 2.0xlO- 3 (^^W j l )(^-) V2 . (61) 
\0.1 pmJ V10" 17 s-V V5 AU/ 

There equations give the radial profiles of (fo, /e, h) for the 
MMSN model at one scale height above the midplane. In ad- 
dition, we need to give the ionization rate ( as a function of 
r. Here, we simply give ( = 10 -17 at r > 3 AU and £ = 10 



at r < 3 AU. The higher value corresponds to ionization by 
cosmic rays and X-rays, while the lower value corresponds to 
ionization by radionuclides. The boundary r = 3 AU is cho- 
sen to approximate the full solution to C(r,z) including these 
ionizing sources (see Figure 2(a) of lO09l) . 

Figure [17] illustrates how the MMSN model is mapped in 
the h-fo plane. This thick solid line in the figure shows the 
radial profiles of fo and h for e = 0.1 and ao = 0.1 pm. The 
line moves upwards in the figure as oq is increased, because 
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fo and h are related as 



f D = 8.4 x 10" 



VO.l fjm) VlO^s" 1 



-6/7 



/j 6/7 (62) 



(this can be directly shown from Equations (l59l and (loTt ) and 
hence /b increases with an for fixed /i. 

Let us see the outcome of fractal dust growth in different 
locations of the disk using the freezeout condition (Equa- 
tion (1571)). Since the condition depends on the three pa- 
rameters (fa, /e, h), the boundary between the growth and 
freezeout regions is a two-dimensional surface in the three- 
dimensional space. However, it will be useful to represent 
the boundary as a single curve in the h-fo plane by relat- 
ing /e to either fa or h. Below, we use the relation = 

l.l(flo/0.1 /im) 11 / 6 /^ 6 obtained from Equations (|59l and 

(go). 

The thick dashed curve in Figure [17] shows below which 
the freezeout condition holds for ao = 0. 1 pro. and e = 0. 1 . For 
this case, we see that the freezeout region covers 1-100 AU 
from the central star. This means that the electrostatic bar- 
rier inhibits fractal dust growth except in an inner region of 
r < 1 AU and an very outer region of r > 100 AU. For com- 
parison, we also show the line (£e / £ Ar) m ax = 3 with the thin 
solid curve (we again use the above relation between /e and 
fa)- This line roughly corresponds to the boundary between 
the growth/freezeout regions predicted by the monodisperse 
theory (see Equation d47|i). Comparing this line with the thick 
dashed curve, we see that the inner region of r < 1 AU would 
be also included in the freezeout region if the bimodal growth 
mode as seen in Section 4 were not absent. From this fact, we 
can expect that the bimodal growth is particularly important 
for dust evolution at small heliocentric distances. It should 
be noted, however, that all these results are dependent on the 
adopted disk model (e.g., laminar disk) a nd param eters (e.g., 
an). We will defer further investigation to Pa per 11 

5.2. Effect of Charge Dispersion 

Up to here, we have assumed that all aggregates with the 
same radius have an equal charge (Q)„. In reality, the charge 
distribution has a nonzero variance, and hence aggregates can 
have a negative charge smaller than the mean value. Here, we 
show that the charge dispersion hardly affects the emergence 
of the total fre ezeou t. 

As shown in O09, the charge distribution for aggregates of 
size a is well approximated by a Gaussian distribution with 
variance (see Equation (24) of O09) 



(SQ 2 



1 + 9 

2 + 9 



ak B T. 



(63) 



In principle, it is possible to fully take this effect into account 
by averaging the collision kernel K over all Q\ and Qi- How- 
ever, the average cannot be written in a simple analytic form. 
For this reason, we simply estimate the effect of the charge 
dispersion as follows. Clearly, the effect of the charge dis- 
persion is significant only if (SQ 2 ) a is much larger than (Q) 2 . 
Using Equations ©, (flBT l. and d63l ). the ratio of (SQ 2 ) a to 
(Q) 2 can be written as 



(SQ 2 



1 + 9 



02)2 2(e E )(2+wy 



(64) 



where (£e) is the electrostatic energy for Q\ = Qi= (Q)a- 
Since 1/2 (1 +*)/(2+ *) sC 1 for all 9, we find that the 
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FIG. 18. — Comparison of the temporal evolution of (N) (left panel) (N) m 
(right panel) between different velocity dispersion models. The thick curves 
are the results for a modified velocity dispersion model (see Section 5.2), 
while the thin curves are the same as the curves showing in the upper panel 
of Figure [14] The thick and thin curves are very similar (indistinguishable 
for h = 10~ 4 5 in the left panel), meaning that the dependence on the velocity 
dispersion is weak. 

ratio (5Q 2 ) a /(Q)l is of an order of (£ E }~ 1 ■ We also find that 
{5Q 2 )a/ (Q)a decreases as dust grows because (£e) increases 
with N. 

Using Equation d64l ). let us consider whether the freeze- 
out criterion (Equation ( f57l i) is affected by the presence of 
the charge dispersion. With the charge dispersion ignored, 
the freezeout criterion is given by (£e)(Nd) > 6. If this con- 
dition holds, we find from Equation d64b that (SQ 2 ) a (N D ) < 
0m(Q)a(N D )] 2 (l + ¥)/(2 + 9) < 0m(Q)a(N D )] 2 . This 
means that the "true" value of £e(Nd) (i.e., the value with 
the charge dispersion taken into account) is not much dif- 
ferent from the "approximate" value (£e)(Nd) as long as 
(£e)(Nd) > 6. Hence, the charge dispersion hardly affects 
the emergence of the total freezeout. 

5.3. Dependence on the Velocity Dispersion 

In this study, we have assumed that the velocity dispersion 
is thermal (see our probability distribution function, Equa- 
tion ((8]l). This assumption neglects any fluctuation in the drift 
acceleration g. This will be reasonable if g is caused by stellar 
gravity (g = fi^z). By contrast, the validity of this approxima- 
tion is unclear if g is driven by turb ulence (g w u„/tn ). For 



example, recent MHD simulations by Carballid o et al. (2008) 
suggest that g may fluctuate by 10% in MRI-driven turbu- 
lence. To check the robustness of our conclusion, we examine 
how the outcome of dust growth depends on the choice of the 
velocity dispersion. 

Here, we consider the cases where the fluctuation in the 
differential drift velocity is as large as the mean value. We 
mimic this situation by replacing k^T in Equation ^ by 
k%T + M m (Amo) 2 , where Auq is the mean relative velocity 
given by Equation ([9). With the modified velocity distribu- 
tion function, we carried out simulations for four sets of pa- 
rameters as in Section 4.2. Figure [18] compare the evolution 
of (N) and (N) m obtained here with that in Section 4.2 (Fig- 
ure [T41) . We find no significant difference between the two 
results. This should be so since the freezeout occurs while 
Brownian motion dominates over the differential drift (i.e., 
k?,T 3> M M (A«£>) 2 ; see Section 4). Detailed inspection shows 
that the average masses grow slightly faster when the disper- 
sion is added to the differential drift, but this is clearly a minor 
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effect. Hence, we conclude that fluctuation in the differential 
drift velocity hardly affects the outcome of the dust growth. 

5.4. Validity of the Fractal Growth Model 

So far, we have assumed that dust grows into porous (frac- 
tal) aggregates. This assumption is true only when the impact 
energy is so low that compaction of aggregates upon collision 
is negligible. Here, we show that the collisional compaction 
is actually negligible when we consider the freezeout of dust 
growth. 

It has been shown by iDominik & Tielensl (119971) that the 
collisional compaction become effective when the impact en- 
ergy exceeds 32^11, where 



w/o charging 



£roll = 37T 7fl fc; 



;6x 10 



-10/ 



7 



100 erg cm 2 ) V 2A 



V2A/V0. 



a 



1/im 



erg (65) 



is the energy needed for a monomer to roll on another 
monomer in contact by 90 degrees. 7 is the surface energy per 
unit area and is estimated as 25 erg cm" 2 for rocky monomers 
and somewhat higher for icy monomers. £ cr ; t is the critical 
rolling displacement for inelastic rolling and is theoretically 
constrained as > 2A ( IDominik & Tielens 19951). 

As seen in the previous section, the total freezeout oc- 
curs when Brownian motion dominates aggregate collision. 
Hence, the relative kinetic energy between frozen aggregates 
is equal to the thermal energy ~ kgT. Assuming T ~ 100 K, 
the thermal energy is ~ 10~ 14 erg, which is many orders of 
magnitude lower than E m \\. Therefore, collisional compaction 
is negligible whenever the total freezeout occurs. 

Of course, the compaction is no longer negligible when 
the electrostatic barrier is overcome since the drift energy in- 
creases with aggregate mass and finally exceeds E m \\. Investi- 
gation of dust growth after the fractal growth stage is beyond 
the scope of this study. 

5.5. On the Role of Porosity Evolution 

As shown in the previous subsection, it is valid to assume 
the fractal dust growth whenever we focus on the freezeout 
of dust growth. However, it has been still unclear whether 
the freezeout occurs even without the porosity evolution. In- 
deed, in previous studies on dust coagulation, it is com- 
mon to ignore the porosity evolution and m odel aggregates 
as spheres of a fixed internal density (e.g., [ Weidenschilling 
1 977[ iNakagawa et al.ll 1 98 ll iTanaka et alj|2005t iBrauer et all 
2008). To fully understand the robustness of the freezeout, we 
will discuss how the growth outcome changes if we adopt the 
compact aggregate model. 

5.5. 1 . Drift and Electrostatic Energies for Compact Dust Particles 

It is straightforward to write down the dimensionless ener- 
gies £0 and £e for the compact model. Since 1Z = N 1 ^ and 
A = N 2 ^ for compact particles, Equations (f2~6b and (|27T i are 
now replaced by 



£d - /d 



NjN 2 
Ni+N 2 



N 2 
A~ 2 



NiN 2 
Ni+N 2 



N, 



and 



! 7^l7^ 2 
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(66) 



(67) 





mass N 

FIG. 19. — Evolution of the mass distribution function J-(N) for compact 
dust models. The gray curves show the snapshots of N 2 J^(N) at T = 10 , 
10 1 5 , . ..,10 3 , 10 3 ' 3 , 10 3 ' 7 (from left to right). The crosses and open cir- 
cles indicate (N) and (N) m at different times, respectively. Parameters 
C/b,/fi,*oc) are set to (10" 5 , 10, 10 05 ). 

respectively. Note that e identically vanishes here by the defi- 
nition of the compact dust model. 

5.5.2. Simulations 

Using Equations d66b and ( t6Tb instead of Equations d26l l 
and d27j >, we have carried out simulations for several sets of 
(fo, fE,h,^>oo). Figure[l9lshows the results for the uncharged 
case (h = 0) and three charged cases (h = 10~ 2 5 , 10~ 4 , 10~ 5 5 ) 
with fixed (//?,/£, *oo) = (10 -5 , 10, 10° 5 ). Note that the val- 
ues of (/dj/ej^oo) are the same as those for the examples 
shown in Sections 4. 1 and 4.2. 

Without charging, the outcome of dust growth is qualita- 
tively similar to that for the porous model (see the upper 
panel of Figure fl2l) . Namely, we see power-law growth at 
early times (T < 10 3 ) and exponential growth at later times 
(T > 10 3 ). One important difference is that the exponential 
growth begins at a lower mass than in the porous case. As 
already mentioned in Section 3, the exponential growth is an 
indication that the differential drift takes over Brownian mo- 
tion in the relative velocity between particles. In the porous 
model, the drift velocity of aggregates increases only slowly 
with mass, because the fractal dimension is close to 2 and 
hence the mass-to-area ratio N/A is nearly insensitive to N. 
In the compact case, by contrast, the drift velocity increases 
with ( Auo oc N/A oc 1//3 ). For this reason, the drift motion 
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FIG. 20. — Outc ome of numerical simulations for compact sphere mod- 
els (see Figure [16] for porous aggregate models). The crosses (x) show the 
parameters for which both (N) and (N) m stop growing at N a Np ("total 
freezeout"), while the open circies (o) indicate where both (N)„, and (JV) 
continue growing ("unimodal growth"). The black dashed curve shows the 
boundary below which £e(No) exceeds 6. For comparison, the boundary for 
porous models (the dashed curve in the left panel of Figure [16) is shown by 
the gray dashed curve. 

takes over Brownian motion (Am oc N~ l l 2 ) at lower N than in 
the porous case. 

The difference mentioned above consequently influences 
the outcome of dust growth with charging charging (the gray 
curves in Figure \\%. We see that the total freezeout does not 
occur at h = 10~ 4 as it does in the porous case. This is because 
of the faster increase in the differential drift velocity men- 
tioned above. In fact, the electrostatic energy also increases 
faster than in the compact case because of the faster decrease 
in the total projected area Aot and capacitance C to t- However, 
this effect is small compared to the faster increase in the ki- 
netic energy. Therefore, we can say that the compact dust 
growth is resistive to the freezeout. Note that the compact 
growth is not free from the occurrence of the freezeout; in 
fact, we observe the freezeout for a higher-/i case, h = 10~ 2 5 . 

We see that the mass distribution for h = 10~ 4 splits into two 
peaks. However, the evolution is qualitatively different from 
what we call bimodal growth in the porous case. The differ- 
ence is that the low-mass peak gets continuously depleted as 
the high-mass peak grows towards higher N. This occurs be- 
cause the high-mass particles acquire arbitrarily high drift ve- 
locities as they grow. For the porous dust model, we have seen 
that the the impact energy for highly unequal-sized collisions, 
£k,\2, is nearly independent of the mass A^i of the heavier par- 
ticle (see Section 4.2). In the compact model, by contrast, the 

impact energy is approximately given by £k\2 ~ 1 + /dNiN^ 3 
(which directly follows from Equation d66l ) with Ny 3> Ay, 
and this increases with A^. However, the electrostatic en- 
ergy £^.12 ~ /bC^/^oo) 2 ^ is independent of Ny as is in the 
porous case. Hence, we find that a high-mass particle with 
sufficiently large Ny can capture smaller particleJJ 

8 Strictly speaking, the decrease in the number of low-mass particles leads 
to the increase in VP (see Section 4.2), and hence proceeds in a way that £ E yi 
balances with £ K yi until 4" reaches ^oo- 



5.5.3. Freezeout Criterion for the Compact Dust Model 

Figure [20] summarizes the results of the simulations for 
compact dust models. The crosses and open circles indicate 
the sets of parameters for which we observe total freezeout 
and unimodal growth, respectively. The gray dashed curve 
shows the boundary below which the freezeout condition sat- 
isfies for the porous model, i.e., the black dashed curve in 
Figurdl6l We see that the compact growth results in the 
freezeout in a more restricted region of the parameter space 
than the porous growth. It is clear that the compact model is 
less conducive to the freezeout compared to the porous model. 

To obtain a freezeout criterion for the compact model, it 
is useful to introduce £q and £e written as a function of 
a single mass N rather than Ny and A^- as done for the 
porous model. There is no difficulty in evaluating £e as- 
suming that the particles are monodisperse, i.e., Ny = Ni = N. 
Using A tot = A/N = AT 1 /} an d c tQt = C/N = AT 2 / 3 , we have 
9 = /i^ooA^. Thus, the electrostatic energy for monodisperse 
compact particles is 



£ E = k[ l+{hN) ^] 



-2.5 . 



N 



1/3 



(68) 



In contrast, we would obtain no meaningful expression for £0 
within the exact monodisperse assumption because £0 identi- 
cally vanishes for N\ = N2= N. For this reason, we will sim- 



ply replace NyNz/(Ny +N2) 
Equation d66l ) to get 



N and |Wj 1/3 -N^ | 



in 



£d = /dN : 



•5/3 



(69) 



As done in Section 3.1, we can define the drift mass No by 

£d(Nd) = 1; using Equation d69l ), we have No = fo^- Hence, 
the critical energy £e(Nd) for the compact model can be ex- 
plicitly given as a function of (fo,fE,h) by 



£e(N d ) ■■ 



f E 



3/5,-0.8 



-25 



f~D 



1/5 



(70) 



Let us examine whether the condition for the freezeout is 
well described by the value of £d(Nd) as is for the porous 
cases. The black dashed curve in Figure [20] shows the line 
where £o(Nq) for the compact model is equal to 6. For com- 
parison, the line £o(Nd) = 6 for the porous case (i.e., the 
dashed curve in the left panel of Figure [16]) is also shown by 
the gray dashed curve. We find that the black line successfully 
explains the boundary between the freezeout and unimodal 
growth regions. Hence, we conclude that the freezeout crite- 
rion for the compact model is again given by Equation ( fSTb if 
only we use Equation < T70b for £o(Nd)- 

To summarize this subsection, we have investigated how the 
growth outcome changes if one adopts a compact dust model. 
We confirmed that the total freezeout does occur even in the 
compact dust growth. This means that a fractal dust model is 
not a prerequisite for the emergence of the freezeout. How- 
ever, this does not mean that the porosity evolution is negligi- 
ble when we analyze the effect of electrostatic barrier against 
dust growth. As shown above, the compact model makes dust 
growth more resistive to the freezeout because the differential 
drift takes over Brownian motion at a lower mass. Therefore, 
the porosity evolution must be properly taken into account in 
order not to overlook the significance of the electrostatic bar- 



rier. 
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In this paper, we have investigated how the charging of dust 
affects its coagulation in weakly ionized protoplanetary disks. 
In particular, we have focused on the effect of the du st siz e 
distribution, which was ignored in the previous work (O09). 
We have us ed the p orous (fractal) aggregate model recently 
proposed by OTS09 to properly take into account the porosity 
evolution of aggregates. 

To clarify the role of size distribution, we have divided our 
analysis into two steps. As the first step, in Section 3, we have 
presented a general analysis on the coagulation of charged ag- 
gregates under the monodisperse growth approximation. The 
monodisperse approximation allows us to define several use- 
ful quantities, such as the maximum energy ratio {£k/ £E)max, 
the drift mass No, and the freezeout mass Nf- We have shown 
that, if the maximum energy ratio {£k / £ £)max is larger than 
unity, the monodisperse growth s talls (or "freezes out") at 
mass N ps Nf, as was predicted by O09. 

As the second step, in Section 4, we have calculated 
dust co agulation using the extended Smoluchowski method 
(OTS09) to examine how the outcome changes when the size 
dispersion is allowed to freely evolve. We find that, under cer- 
tain conditions, the electrostatic repulsion leads to bimodal 
growth, rather than total freezeout. This bimodal growth is 
characterized by a large number of "frozen" aggregates and 
a small number of "unfrozen" aggregates, the former control- 
ling the charge state of the system and the latter growing larger 
and larger carrying the major part of the system mass. 

Based on the results of our numerical simulations, we have 
obtained a set of simple criteria that allows us to predict how 
the size distribution evolves for given conditions (Section 4.3; 
Tabled. These read: 

• If £e(Nd) > 6, all aggregates stops growing before the 
systematic drift dominates their relative velocities (total 
freezeout). The outcome is a nearly monodisperse dis- 
tribution of frozen aggregates with typical mass s» Nf- 

• If £e(Nd) < 6 and < *oo/4, a large number of ag- 
gregates stop growing, but the major part of dust mass 
within the system is carried by a small number of ever- 
growing aggregates {bimodal growth). 

• If £ E (N D ) < 6 and > *oo/4, all aggregates con- 
tinue growing with a single-peaked size distribution 
(unimodal growth). 

The second case includes situations where aggregates cannot 
continue growing in the monodisperse growth model. Thus, 
the size distribution is an important ingredient for the growth 
of dust aggregates beyond the electrostatic barrier. 

We emphasize again that our analysis assumed fractal evo- 
lution of dust aggregates. This assumption is valid only when 
the collision energy is so small that collisional compactio n 
is negligible (Dominik& Tielens 1997; Suvamaet al. 2008). 
We have proven that the collisional compaction is indeed neg- 
ligible as long as the total freezeout is concerned since the 
freezeout always occurs when Brownian motion dominates 
aggregate collision (Section 5.4). It should be noted that most 
theor e tical studies on dust coagulation (e.g., Nakagawa et al. 
11981b iTanaka et al] 12005b iBrauer et al.1 120081) have ignored 
the porosity evolution and modeled aggregates as compact 
spheres. However, we have found that such simplification 
leads to underestimation of the electrostatic barrier because 
compact spheres are frictionally less coupled to the gas and 
hence have higher drift velocities than porous aggregates of 
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FIG. 21 . — Mass-to-area-ratio B = N/A versus monomer number N for 
numerically created BCCA clusters. The thin solid curves show 20 sam- 
ples, while the thick solid curve indicates the ave rage over 100 samples. The 
dashed curve shows Minato's formula (Equation \2\Y ). 

the same mass (Section 5.5). Therefore, the porosity evolu- 
tion must be properly taken into account when considering 
the electrostatic barrier against dust growth in protoplanetary 
disks. 

In iPaper PL we apply our growth criteria to particular pro- 
toplanetary disk models to investigate the effect of the elec- 
trostatic barrier in the early stage of planet formation. 

The authors thank the anonymous referee for the many 
comments that greatly helped improve the manuscript. S.O. is 
supported by Grants-in-Aid for JSPS Fellows (22 ■ 7006) from 
MEXT of Japan. 

APPENDIX 
NUMERICAL ESTIMATION OF THE AREA 
DISPERSION 

Let us consider two groups of porous aggregates each of 
which is characterized by aggregate mass Nj(j = 1,2). In ei- 
ther group, aggregates have different values of the projected 
area Aj. Therefore, the projected area, or the mass-to-area 
ratio Bj = Nj/Aj, of an aggregate randomly chosen from the 
j-th group can be regarded as a stochastic variable. The aver- 
age of the quantity \B\-B2\ 2 over all possible pairs is given 
by 

|Bi-B 2 | 2 = ^{N^-BiNit+Yj&Wi^ 

7=1,2 

= \B(N l )-B(N2)\ 2 + <Nj) 2 B(Nj)\ (71) 

7=1,2 

where B(Nj) and 8B 2 (Nj) are the statistical average and vari- 
ance of B for aggregates of the y'-th group, and e(A0 = 

SB 2 (N) 1/2 /B(N). Note that we have assumed that B { and B 2 
are uncorrected, i.e., = B(N{)B(N2)- Equation ( fTTT i re- 
duces to Equation (l23l if e(A0 is independent of N. In this 
appendix, we estimate e(N) using numerically created BCCA 
clusters. 

We have performed 100 BCCA simulations and obtained 
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F I G 22. — Normalized area dispersion e = SB 2 fB for sample BCCA 
clusters. The solid and dashed curves are obtained by averaging 100 and 50 
samples, respectively. 



the relation between A and N for each run. Since the pro- 
jected area of an aggregate generally depends on the choice 
of the projection angle, we determined it as the average over 



15 randomly chosen orientations. Figure |2TI shows the mass- 
to-area ratio B versus monomer number N for 20 samples as 
w ell as the average B over 100 samples. The area formula 
of Mi nato et al.l (2006). Equation ( |2T1 i. is also plotted to show 
that B is consistent with the finding of Min ato et alj ((2006). 

Figure [22] shows the ratio e(A0 obtained from 100 samples. 
For 10 < N < 10 6 , e(A0 is of an order of 10 _1 and increases 
very slowly withAf. Therefore, e(A0 can be well approximated 
as a constant 10 _1 . To check the convergence, we compute 
e(N) using 50 of the samples. The small difference between 
the two curves means that the statistical error due to the finite 
number of samples is negligible. 

Figure|22]implies that e(N) may be considerably larger than 
10" 1 for N 3> 10 6 . However, it should be noted that the above 
clusters has been formed through collisions between identi- 
cal clusters. In reality, an aggregate in an ensemble collides 
with aggregates of various sizes. The most probable are col- 
lisions between aggregates of very different B, since the col- 
lision probability is proportional to |Bi— Bz\. This effect gen- 
erally cause the decrease in SB 2 , and hence the decrease in e. 
Therefore, the value of e estimated here should be regarded as 
the upper limit of the actual values. 
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